This invention relates to methods for generating and processing simultaneous vibratory seismic data.
In seismic prospecting using simultaneous vibratory techniques, a series of seismic energy sources or vibrators are employed to transmit seismic signals into the earth. Part of these signals are reflected from interfaces between subterranean strata, and/or refracted within strata, back to the surface of the earth, where they are detected by one or more receivers. The time taken for a signal to pass from a particular vibrator to a particular receiver gives an indication of the length of travel of the signal between that vibrator and that receiver, from which the structure of geological formations may be deduced.
Vibrators generally consist of a base plate in contact with the ground, with a large weight applied to this plate. The seismic energy is transmitted to the ground by applying a vibratory force to the plate. Thus if a constant weight of 30 tonnes is applied to the metal plate, a vibration with amplitude 20 tonnes can be applied to the vibrator, ensuring that there is always a positive force against the ground.
Many modern vibratory seismic surveys are performed by simultaneously imparting energy into the earth from multiple source locations, so that each receiver will detect refracted and reflected energy which has been emitted by the whole series of vibrators. The data recorded at each receiver must then be processed so that the signal due to each individual vibrator can be separated out.
This is generally achieved by each vibrator performing multiple xe2x80x9csweepsxe2x80x9d or xe2x80x9cshotsxe2x80x9d, where the phase of the signals emitted by the vibrators are varied between vibrators and between shots. In its simplest form, this can be illustrated by the case of two vibrators, twice operated simultaneously. If they are operated in phase with each other for the first sweep, but 180xc2x0 out of phase for the second sweep, the receiver will record two signals. These signals may be added together to determine the signal arriving from the first vibrator, or subtracted to determine the signal arriving from the second vibrator.
In general, in order for the data arriving at each receiver to be decoded so that the contribution of each vibrator can be determined, there must be at least as many shots as there are vibrators. For each shot by each vibrator, a waveform is applied to that vibrator by the control mechanism. This waveform is normally a constant amplitude swept-frequency xe2x80x9cchirpxe2x80x9d, with tapered ends due to the fact that the amplitude has to be ramped up at the start and ramped down at the end, and is known as the xe2x80x9cpilot sweepxe2x80x9d. In modern seismic operations, the pilot sweep almost always begins at low frequency and finishes at high frequency, and the frequency normally increases linearly with respect to time.
In practice, the waveform actually applied to the ground by the vibrator is not quite the same as the waveform applied to the vibrator. Deviations from the pilot sweep are inevitable. These deviations have two characteristic forms:
1. As well as applying a force at the desired frequency (known as the fundamental), the vibrator also applies a force at integer multiples of that frequency (known as harmonics).
2. The force applied at the fundamental frequency deviates in amplitude and/or phase from the pilot sweep.
At each receiver, in order to separate out the signal from each vibrator, it is necessary to use some approximation of the signal provided by each vibrator. In the past, two methods have chiefly been used.
The first assumes that the vibrator force follows the pilot sweep exactly, and separates out the earth response due to the vibrator at each vibrator location accordingly. If such a method is used, the data at each receiver, or geophone, can be correlated with a single representative pilot sweep before separation. No measurements need to be made on the vibrator, and the procedure is robust. Non-linear effects such as a non-linear earth response or data clipping during acquisition will not leave significant artifacts.
The main disadvantage with this method is that deviations of the force applied at the fundamental from the pilot sweep can lead to significant cross-coupling of the vibrators, i.e. the signal from one vibrator being ascribed to another. This occurs if each vibrator has a different deviation from the pilot sweep. This method can also lead to timing errors for the individual vibrators.
The second known method uses a measurement of the force applied to the earth by the vibrator. Thus, instead of using the pilot sweep, the real force applied by the vibrator is used to derive an inversion operator. This removes the main disadvantage of using the pilot sweep, but itself has other disadvantages. The whole inversion procedure must be performed on uncorrelated data, the entire waveform of each shot must be acquired for each vibrator and each shot, and the method is sensitive to the non-linear effects described above.
In general, the data acquisition and much of the processing are performed separately. The data recorded by the receivers may be partially processed as it is received, but normally either raw or partially processed data is saved onto tapes which are then transferred to a central data processing unit.
According to a first aspect, the present invention provides method of processing seismic data, said seismic data having been obtained by:
performing a plurality of sweeps, wherein each sweep comprises generating seismic signals in the earth using a plurality of vibrators by applying a pilot sweep waveform to each vibrator, each pilot sweep being a waveform of changing frequency;
measuring the force applied to the earth by each vibrator to determine a measured force waveform; and
measuring the seismic signals at one or more locations remote from the vibrators;
said method comprising:
filtering the measured force waveform to remove harmonics of the pilot sweep and thus determining a filtered force waveform;
generating an inversion operator from the filtered force waveform for each vibrator; and
applying said inversion operator to the measured seismic signals to determine the contribution of each vibrator to the seismic signals.
The force waveform applied to the earth may be measured by special sensors, although preferably it is measured by derivation from the signal used by the vibrator controller to control the vibrator output. This may be a weighted sum measurement, consisting of a linear combination of signals from accelerometers on the vibrator.
The filtered force waveform may be determined by the application of a time varying notch filter to the measured force waveform, but more preferably the measured force waveform applied to the ground is filtered by cross-correlation with the pilot sweep followed by application of a time window to remove the harmonics of the pilot sweep. Each measured waveform may preferably be cross-correlated with the same pilot sweep. The time window is preferably applied around zero time.
The measured seismic signal is preferably also filtered by cross-correlation with the pilot sweep followed by application of a time window before the application of the inversion operator.
The pilot sweep preferably comprises a waveform having a substantially constant amplitude envelope with tapered ends, within which the frequency increases with time.
The number of sweeps may be the same as the number of vibrators. In another embodiment, the number of sweeps is greater than the number of vibrators, and the noise on each measured seismic signal is estimated and used in the determination of the inversion operator.
In a preferred embodiment, the invention provides a method of performing a seismic survey, comprising:
performing a plurality of sweeps, wherein each sweep comprises generating seismic signals in the earth using a plurality of vibrators by applying a pilot sweep waveform to each vibrator, each pilot sweep being a waveform of changing frequency;
measuring the force applied to the earth by each vibrator to determine a measured force waveform;
measuring the seismic signals at one or more locations remote from the vibrators; and
processing seismic data comprising the measured seismic signals and the measured force waveform using any of the methods described above.
At least in preferred embodiments, the data can be correlated before inversion, and the cross-coupling and timing error problems are removed. The method is not sensitive to non-linear effects. Further, a greatly reduced quantity of data need be recorded for each vibrator.
If there are the same number of vibrators as sweeps the inversion procedure is straightforward. However, the method may comprise taking more sweeps than there are vibrators.
According to a second aspect, the present invention provides a method of processing seismic data, said seismic data having been obtained by:
performing n sweeps, each sweep comprising generating seismic signals in the earth using nm vibrators in contact with the earth, where 12 is greater than m;
measuring the seismic signals at one or more locations remote from the vibrators to determine the seismic data; and
estimating the noise N associated with the seismic signals for each sweep;
the method comprising:
generating a noise correlation matrix P from the noise N associated with each sweep;
generating a nxc3x97m matrix M from the signal provided by each of the is sweeps at each of the m vibrators, each element of the matrix M corresponding to a source signature from one vibrator for one sweep: and
generating an inversion operator Mxe2x8axa5 for application to said seismic data, said inversion operator being determined by
Mxe2x8axa5=(M*P(xe2x88x921)M)(xe2x88x921)M*P(xe2x88x921).
The method preferably also comprises selecting rows from M in order to generate all possible square mxc3x97m sub-matrices of M;
wherein the inversion operator Mxe2x8axa5 is generated by performing a relative weighted sum of the inverses of the square mxc3x97m sub-matrices, where the relative weight applied to each inverse is the product of the elements of P corresponding to the rows omitted from M to generate the sub-matrix, multiplied by the absolute value of the determinant squared of the sub-matrix.
P may be determined by generating the component Pij of P in the ith column and jth row using Pij= less than NiNj* greater than .
Preferably, the inversion operator Mxe2x8axa5 is determined by       M    ⊥    =                              (                                    M              *                        ⁢                          P                              (                                  -                  1                                )                                      ⁢            M                    )                          (                      -            1                    )                    ⁢              M        *            ⁢              P                  (                      -            1                    )                      =                            ∑                      ij            ⁢                          xe2x80x83                        ⁢            …                          ⁢                              (                                          ∏                                  ij                  ⁢                                      xe2x80x83                                    ⁢                  …                                            ⁢                              P                k                                      )                    ⁢                      (                                          "LeftBracketingBar"                                  det                  ⁡                                      (                                          M                                              ij                        ⁢                                                  xe2x80x83                                                ⁢                        …                                                              )                                                  "RightBracketingBar"                            2                        )                    ⁢                      M                          ij              ⁢                              xe2x80x83                            ⁢              …                                      (                              -                1                            )                                                            ∑                      ij            ⁢                          xe2x80x83                        ⁢            …                          ⁢                              (                                          ∏                                  ij                  ⁢                                      xe2x80x83                                    ⁢                  …                                            ⁢                              P                i                                      )                    ⁢                      (                                          "LeftBracketingBar"                                  det                  ⁡                                      (                                          M                                              ij                        ⁢                                                  xe2x80x83                                                ⁢                        …                                                              )                                                  "RightBracketingBar"                            2                        )                              
Where Mij . . . is the mxc3x97m sub-matrix.
There may be one more sweep than vibrator, in which case the inversion operator Mxe2x8axa5 may be determined by       M    ⊥    =                    ∑                  i          =          1                n            ⁢                        P          i                ⁢                              "LeftBracketingBar"                          det              ⁡                              (                                  M                  i                                )                                      "RightBracketingBar"                    2                ⁢                  M          i                      (                          -              1                        )                                              ∑                  i          =          1                n            ⁢                        P          i                ⁢                              "LeftBracketingBar"                          det              ⁡                              (                                  M                  i                                )                                      "RightBracketingBar"                    2                    
where Mi is the mxc3x97m sub-matrix omitting the ith row of M.
The noise power estimates may preferably be made using quadrature mirror filtering, although alternative methods such as discrete cosine transforms or discrete Fourier transforms may also be used.
In a preferred embodiment, the invention provides a method of performing a seismic survey, comprising:
performing n sweeps, wherein each sweep comprises generating seismic signals in the earth using m vibrators in contact with the earth, where n is greater than m;
measuring the seismic signals at one or more locations remote from the vibrators to determine seismic data;
estimating the noise N associated with the seismic signals for each sweep; and
processing the seismic data using any method described above.
According to a third aspect, the present invention provides a method of generating seismic data, comprising:
performing a plurality of sweeps, each sweep comprising generating seismic energy using a plurality of vibrators, each vibrator operating at a different phase, the power applied to each vibrator for each sweep being substantially the same;
wherein the phases are cycled through the vibrators, in such a way that for n vibrators labelled 1 to n, for each sweep after the first, vibrators number 2 to n use the phases of the vibrators 1 to nxe2x88x921 respectively from the previous sweep:
and wherein if a matrix is determined from the phases of the vibrators, with each column corresponding to one sweep and each row corresponding to one vibrator, the magnitude of the determinant of the matrix is the same as the dimension of the matrix if the phases of the vibrators are written as complex numbers.
The number of sweeps and number of vibrators may be the same. Alternatively, there may be more sweeps than vibrators, in which case a matrix for the phases of the vibrators may be determined as if there were as many vibrators as there are sweeps, and columns deleted from said matrix so that there are as many columns as there are vibrators.
Preferably, there is one more sweep than there is vibrator.
This method of determining the phases of the vibrators minimises the ambient noise at each receiver by keeping the size of the inversion operator as small as possible.
In preferred embodiments, a seismic survey may be conducted using seismic signals generated using the methods just described and processing data using any of the methods described above.
According to a fourth aspect, the present invention provides a computer readable medium carrying a computer program arranged to perform any of the methods described above.